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The strongly nonlinear regime of Dyakonov-Shur insta- 
bility is studied analytically and by computer simulations. 
The instability leads to high-amplitude stationary oscillations 
caused by shock wave formation in the FET channel. The nu- 
merically found shape and speed of the shock waves are well 
described in terms of the conventional hydrodynamics. The 
amplitude of stationary oscillations is determined. 

PACS numbers: 72.30-f-q, 73.20Mf 



I. INTRODUCTION 

Dyakonov and Shur flhave described a new mecha- 
nism of plasma wave generation in an ultrashort-channel 
(ballistic) Field Effect Transistor. It was shown that 2D- 
electrons in a BFET can be described by equations anal- 
ogous to the hydrodynamic equations for shallow wa- 
ter, the plasma waves being similar to shallow water 
waves. Therefore, various hydrodynamic phenomena ob- 
served in shallow water may occur in the 2D-electron 
fluid. In Ref. it was demonstrated that the non- 
symmetric boundary conditions with fixed source poten- 
tial and drain current lead to instability of the stationary 
state with dc current. The FET channel is a resonator 
cavity, with plasma wave generation mechanism similar 
to that known for self-excitation of jets and organ pipes. 
Both the eigenfrequencies and the instability increment 
were found for dissipationless electrons. The plasma wave 
generation in short channel FETs is very important for 
using these devices as high-power sources in th^ terahertz 
frequency range. The experimental results o' ™ found 
for the response of GaAs HEMT inj-jthe detector mode 
agrees with theoretical predictions, o Then, the recent 
experiments concerning GHz plasma waves in 2DEG 
make an observation of " shallow water instability" to be 
more feasible. 

In Ref. , nonlinear evolution of the instability was 
considered with account of both the electron fluid viscos- 
ity related to electron-electron scattering and the friction 
due to electron scattering by phonons and(or) impuri- 
ties. The instability has a threshold found for the case 
of small viscosity and friction. B As expected, the insta- 
bility leads, instead of chaotic motion, to the establish- 
ment of low-amplitude stationary oscillations, provided 
that the threshold current is exceeded by a small enough 
margin. In accordance with the theow of weak turbu- 
lence in conventional hydrodynamics , the amplitude 
is proportional to the square root of the current near 



the instability threshold. As predicted in Ref. at 
high enough increments the instability may result in the 
formation of stepJike electron concentrations and veloc- 
ity distributions B similar to shock waves in the conven- 
tional hydrodynamics. In Ref. |^ the shock wave forma- 
tion for non- viscous flow in the FET channel for the case 
of symmetric boundary conditions was investigated both 
analytically and numerically. 

In Ref. [|lOj, the instability increment and threshold 
were calculated in the special case of strong friction and 
zero viscosity. The electron scattering results in the nar- 
rowing of an instability region, with the instability in- 
crement decreasing. The instability vanishes at a certain 
critical value of friction. The analytical results are in 
excellent agreement with computer simulations. 

Neglecting the friction of the electron fluid, we study 
the strongly nonlinear instability mode both analytically 
and using computer simulations. First, we check the reli- 
ability of our numerical method in the analytically solv- 
able case of a linear instability mode. Then, we demon- 
strate that the instability leads to high-amplitude sta- 
tionary oscillations caused by shock wave formation in 
the FET channel. The shape and speed of the numer- 
ically found shock waves are well described in terms of 
the conventional hydrodynamics. Finally, the amplitude 
of stationary oscillations is determined. 



II. ANALYTICAL APPROACH 

In Ref. |l| , it was shown that 2D electrons in a BFET 
can be described by the equations analogous to the hy- 
drodynamic equations for shallow water: 
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Here, the voltage swing, U = Ugc — Ut, corresponds to 
shallow water level, Ugc is the local gate-to-channel volt- 
age, and Ut is the threshold voltage. Then, V is the 
local electron flux velocity, m is the electron effective 
mass, K is the viscosity related to electron-electron scat- 
tering. Eq. (|l|) is similar to the continuity equation in 
hydrodynamics, taking into axjcount the gradual channel 
approximation relationship O 



CU 

e 



(3) 



1 



Here, Ug is the surface electron concentration, C is the 
capacitance per unit area. Then, Eq.(|^) is analogous to 
the Euler's equation in which we neglect the pressure gra- 
dient term. B Note, Eqs.(|l[]|) represent, in fact, the sim- 
plified set of equations. In general, one have to solve the 
exact Poisson's and Navier-Stocks hydrodynamic equa- 
tions. Actually, within gradual channel approximation 
the higher order terms in the plasmon dispertion law O 
as well as the edge-effects contribition to the electric field 
are neglected. 

In Ref. 0, it was demonstrated that the steady elec- 
tron flow is unstable under the following boundary con- 
ditions: 



C/(0,t) = C/o, U{l,t)V{l,t) = J/C, 
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where Uq is the fixed source potentials; = 0), and J 
is the current per unit gate width, fixed at the drain 
(x = I). In the steady state, V,U are constants: U = 
Uo,V = Vq = j/(CUo)- Neglecting the viscosity, we can 
easily solve Eqs.(^[^) linearized with respect to small per- 
turbations of the steady state flow (see Appendix). As 
shown in Ref. [Q, the instability results in plasma wave 
generation, with the FET channel playing the role of a 
resonator cavity. Both the eigen frequencies and the in- 
stability increment (see Eq.A-II) depend on M = Vq/S, 
where M is the analog of the Mach number in hydrody- 
namics and, S — {eUo/mY^'^ is the plasma wave velocity. 
The instability exists at < M < 1. 

We now investigate the strongly nonlinear instability 
mode, when shock waves are formed in the FET chan- 
nel. As well known in the conventional hydrodynamics, 
the mass and momentum conservation laws determine 
the shape and variables on both sides of the shock wave. 
With dimensionless coordinate rj = x/l, time r — tS/l, 
and variables u = U/Uq,v = V/S introduced, Eqs.(0||) 
yield 



du 
97 



d{vu) d 



dr 



drj 



d{vu) 
df] 



= 0, 



(5) 
(6) 



where k = x/ISUq is the dimensionless viscosity. 

We can easily find a self-consistent solution of Eqs. (^|j|) 
which corresponds to a shock wave propagating in an infi- 
nite medium at fixed velocity, c, of the wave front. Simple 
integration of Eqs.(^,||) over the interval containing the 
shock wave discontinuity yields 
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vu{v — c) H = K— 
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where j, q are the current and momentum fluxes which 

are continuous in the shock wave. With Eq.(0), the shape 
of the shock wave can be easily derived as follows: 
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where Vr = v' /j^^"^, v' — v — c is the relative velocity of 
the electron flow with respect to the shock wave front, 
const = {jc — q)/j^^^- It is noteworthy that for a shock 
wave propagating in the inflnite medium the dimension- 
less coordinate, ^ = rjj / k = xJSUo/k, does not contain 
the sample length scale I. Then, according to Eq.(||) the 
spatial distribution of the shock wave velocity and poten- 
tial, u(^) = j^^^ /vriS,)^ universal functions of ^. For 
a low-amplitude shock wave, solving Eq.(P) yields 



U1+U2 , U2-U1 
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Here, the indices 1 and 2 label, respectively, the values 
of the hydrodynamic variables in front of, and behind 
the shock wave (see Fig.^,b), in accordance to Ref. [Q. 
Then, S = „ -^tT"^ — r is the front width of the low- 
amplitude shock wave. The smaller the amplitude, the 
wider the front. 

We emphasize that Eq.(^ provides that the energy 
conservation law is already satisfied. The energy dissi- 
pated within the discontinuity is equal to the difference 
of energy flux densities j{u+ (v')^ /2) on both sides of the 
shock wave. Moreover, Eq.(|^) determines the direction in 
which the variables in the shock wave change as U2 > ui, 
then gives the jump conditions known in the conventional 
hydrodynamics. \3 Indeed, for a shock wave in the infi- 



nite medium ^ 
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— 0, and, therefore the components 



Ml. 2 and vi,2 are related by 
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where 2 = ^'i,2~c is the relative velocity of the electron 
flow in front of, and behind the shock wave, respectively. 

We now seek the velocities and the potential on both 
sides of the shock wave and the speed of the shock wave 
front, taking into account Eq. (p^) and the dimension- 
less boundary conditions ui = l,M2i'2 = M specified 
by Eq.(^). In this case, however, the number of alge- 
braic equations is less than the total number of variables. 
Consequently, we set the relative shock wave amplitude, 
U21 = U2/U1, as a free variable. Finally, the speed of the 
shock wave front, c± , and the electron flow velocity at 
the source, vf , are given by 



M 

U21 
M 



'Ji = ±v{l- U21). 

U21 



(11) 

(12) 



Here, the signs correspond to shock wave propagating in 
the channel upstream (+) and downstream (-), respec- 
tively. Then, v = ^ ^2^"^ ) shock wave front 
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speed for a stagnant electron fluid (i.e. that with M = 0). 
It should be noted that the speed of the low-amplitude 
( U21 — 1) shock wave front coincides with the plasma 
wave propagation velocities Mi 1, in agreement with the 
conventional hydrodynamics. It will be recalled that 
Eqs.(^|l]) are vahd when U21 > 1- 



III. COMPUTER SIMULATIONS 

Let us first estimate the order of magnitude of the di- 
mensionless viscosity k. In a highly non-ideal 2DEG, 
with the thermal, Bohr, and Fermi energies of the same 
order of magnitude, the viscosity of the electron fluid is 
given Gf by VpXee ^ h/m, where Vp is the Fermi ve- 
locity. Age ns is the mean free path for electron- 
electron collisions. For Uq ^ 0.5V, Us ^ lO^^cm"^, 
and submicrometer( 1 — 0.2 /im ) AlGaAs/InGaAs based 
FET at T = 77K we obtain k = XeeVp/lS « 0.01. Ac- 
cordingly, we further assume that k <C 1. 

We now compare our analytical results with computer 
simulation data. The numerical-method used is the well- 
known Broilovskaya algorithm 113 adapted for a system 
with two boundaries. In Ref. a detailed description 
of the numerical procedure based on Eqs.(^j6|) was given. 
As was underlined in Ref. , a finite-difference approx- 
imation methods always involve a numerical error which 
could be associated with a diffusion-like term similar to 
that in Eq.(||). Consequently, we further assume that 
within our numerical method there exists a nonzero par- 
asitic "numerical viscosity" Kn- We will demonstrate that 

Kn < 1. 

In the present paper we use two different numerical 
approaches. Initially (case (i)), we omit viscosity term 
(k = 0) in Eqs.(||,^, and then solve them under the 
boundary conditions specified by Eq. (^) . We use the 
first mode of plasma wave oscillations ( Eqs.(A-III) at 
A = 10-4,^ 



= l,r = 0) as the initial condition. In 
Fig.|l|, the computer simulation data are presented as the 
time dependence of the velocity at the source, w(0, r), in 
relation to the drain potential u(l,r) for different Mach 
numbers. In the same figure, we plot the solution given 
by Eq.(A-II]), found within the linear approximation. 
The computer simulation data are well described by the 
analytical theory in the linear regime of the instability 
evolution for < M < 1. 

We now consider the strongly nonlinear instability 
regime. Instead of a chaotic motion, as expected, the 
instability leads to stationary oscillations irrespective of 
the initial conditions. The oscillation amplitude depends 
on the Mach number. As an example, in Fig.^ we plot 
four different phases clearly defined in stationary oscilla- 
tions. In the first step (Fig.^,a), the shock wave appears 
at the drain and moves therefrom toward the source. Af- 
ter subsequent reflection from the source the jump propa- 
gates upstream the channel and then disperses (Fig.^,b). 
The dispersion of the shock wave front is a generic phe- 



nomenon known in the conventional hydrodynamics. 
In other phases, there are relatively smooth distributions 
moving downstream and upstream the channel. It is note- 
worthy that the shallow rapid oscillations near the dis- 
continuity( see Fig.^,a ) are numerical artifacts that al- 
ways appear in finite-diflference approximation methods 
for a non- viscous flow. 113 

Let us now examine in more detail the shock wave 
clearly defined (bold line) in Fig.|^,a. It can be seen 
that a discontinuity approaching the source can be de- 
scribed analytically in terms of a shock wave in the infi- 
nite medium. To confirm this, we plot in Fig. p] th e an- 
alytical dependences vf{M,U2i) given by Eq.(|12[). We 
ascertained that the upstream shock wave at the source 
is well described by the theory at 0.05 < M < 1. It 
should be noted that the above description fails to ac- 
count for the dependence of the amplitude of station- 
ary oscillations on the Mach number. Nevertheless, we 
use the computer simulation data and find the ampli- 
tude of stationary oscillations, ufJM), associated with 
the downstream shock wave (Fig.|l|, point A). This de- 
pendence (see Fig.^,a) correlates with that furnished by 
Eq.(A-II) for the instability increment. The higher the 
increment, the stronger the amplitude of stationary os- 
cillations. Using the dependence U2i{M) and Eq. (pi]) , 
we can readily calculate the speed of the upstream and 
downstream shock waves (see Fig.^j) . It can be seen that 
at M ^ 0; 1 the speed of the low-amplitude shock waves 
coincides with the plasma wave propagation velocities 
M±l. 

Up to now, we have used the approach k = 0, with 
the parasitic "numerical viscosity" assumed to be small 
Kn <C 1. Let us examine the opposite situation (case(ii)) 
K„ ^ K 1 with viscous term included in Eq. (|6|) . If this 
case, Eqs.(^,^) has higher order and, consequently, re- 
quires an additional boundary condition other than those 
given by Eq.(U). Fortunately, the extra condition can 
be readily found with the help of the continuity equa- 
tion. Indeed, for a fixed voltage at the source Eq.(||) 

= 0. Previously, this condition was derived 

on the basis of a more complicated reasoning. 
With the above set of boundary conditions, we examined 
numerically the instability evolution of a viscous flow at 
0.005 < K < 0.05 and M < 0.15. As expected, in the lin- 
ear regime of instability, the numerical data coincide with 
those found above for the non- viscous case (i). Then, in 
the nonlinear regime, the viscosity enhancement results 
in a minor damping of stationary oscillations. The am- 
plitude of these oscillations remains nearly constant in 
the actual range of k, and, hence, we can use the simple 
non-viscous approach (i) instead of (ii) in this case. 

To conclude, we estimate the parasitic "numerical vis- 
cosity." In Fig.^, we plot the spatial distributions for up- 
stream shock waves(see Fig|^,a, point A) at fixed Mach 
number M = 0.1 and different viscosities, it will be re- 
called that, for a shock wave in the infinite medium, the 
spatial distribution u(^) scales with ^ ^ k. Consequently, 



yields |i 



Ref. § 
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we can represent the same data re-plotted, e.g. the case 
K = 0.01. The plots are collapsed(see Fig.Q,b). In addi- 
tion, we plot in the same figure the distribution of the 
potential in the shock wave, found from Eq.^. It can be 
seen that the shape of the numerically found shock wave 
is well described by the theory. Using the collapse for- 
malism we can easily estimate the parasitic "numerical 
viscosity" to be k„ = 0.0013. This value is consistent 
with estimations done in Ref. jl^. To the best of our 
knowledge, the exact value of the "numerical viscosity" 
within the Broilovskaya algorithm is unknown. 

In summary, in the nonlinear regime the current in- 
stability leads to high-amplitude stationary oscillations 
caused by shock waves formation in the FET channel. 
The amplitude of stationary oscillations is determined 
using computer simulation data. The numerically found 
shape and speed of the shock waves are well described in 
terms of the conventional hydrodynamics. The parasitic 
"numerical viscosity" is estimated. 
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V. APPENDIX 

Following Ref. Q, we seek the instability by analyz- 
ing the temporal behavior of small perturbations (5V, SU 
superimposed on a steady flow, i.e. V{x, t) = Vq + 
SV{x) ex-p{—iujt), U{x,t) — Uo + SU{x) exp{—iujt). Us- 
ing dimensionless variables introduced in the text and 
non-viscous(K = 0) Eqs.(|l|-H) linearized with respect to 
SV, 5U we obtain 

V = M + ARe [exp(i(fci?7 - f^r)) exp{i{k2r] - nr))] (A-I) 
u = 1 + ARe [exp{i{kiri — ilr)) — exp(i(A:2?7 — ^t))] 

where A is the arbitrary constant dependent on the initial 
conditions, and ~ ujI / S is the dimensionless frequency. 
Then, /ci,2 = ±il/(l ± M) denote the wave vectors for 
the upstream and downstream plasma waves respectively. 
According to Ref. ||l|, both the real and imaginary parts 
of the complex frequency = Q' + iQ" of plasma wave 
generation are given by 



n' = 



(1 



M2) „ (1 - M2) 

Un, n" = ^ L In 



1 + M 



1 - M 



(A-II) 



where n is a n od d integer for \M\ < 1. 

From Eqs.( A-]) we can readily find the time-dependent 
velocity at the source w(0,t) and drain potential u(l,r) 
as follows 



■i;(0, T)=M + 2ARe [exp(-if)r)] 
«(1,.)^1 + 2aL_^ 



(A-III) 

- Re [exp(i(Av3 - fir))] 



where A(p = 7rn(l — M)/2 is the phase shift. We empha- 
size that Eqs.( A-III) describe a helix (see Figj^). 
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FIG. 1. Instability evolution diagram for (a) Al — 0.1, (b) 
M = 0.3, (c) M = 0.5, and (d) M = 0.7. The thinUne repre- 
sents the linear regime of instability, given by Eq.( |A-II]| ). The 
dashed (dotted) asymptotes are given by Eq.(|l^) at U21 > 1. 
The numbers denote the time slice numbers(see Fig.^. 



FIG. 2. Spatial distribution u{£^) within the stationary os- 
cillation period at M = 0.1. 



FIG. 3. Stationary amplitude and the speed c± of the 
shock waves at the source ( Fig |l|, point A) vs the Mach num- 
ber. 
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FIG. 4. The graphic solution of Eq.(H): (a) the potential 
u(^) and velocity Vr{^) distributions in the shock wave (b) the 
plot dvr/d^ vs Vr- Computer simulation data for: (c) spatial 
distribution u{ri) at M = 0.1 and k = 0.005; 0.01; 0.02. 
and (d) the same curves collapsed with respect to k = 0.01. 
Dashed line: analytical result given by Eq.(H) 
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